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We demonstrate the practical feasibility of calculating transport coefficients such as the viscosity of 
liquids completely from first principles using the Green-Kubo relations. Results presented for liquid 
aluminum are shown to have a statistical error of only ca. 5%. The importance of such calculations 
is illustrated by results for a liquid iron-sulfur alloy under Earth's core conditions, which indicate 
that the viscosity of the liquid outer core is not substantially higher than that of typical liquid 
metals under ambient conditions. 



PACS numbers; 66.20.+d, 7f.f5.Pd 

Since the pioneering work of Car and Parrinello |^] 
first principles molecular dynamics (FPMD) has become 
a widely used technique for investigating condensed mat- 
ter. For liquids and solids in thermal equilibrium many 
quantities of interest are readily calculated, including 
thermodynamic functions, structure factors, radial dis- 
tribution functions, diffusion coefficients, bond lifetimes, 
etc. But transport coefficients such as viscosity and ther- 
mal conductivity are computationally more demanding, 
and have been little studied by FPMD. We demonstrate 
here the feasibility of calculating such quantities with 
useful accuracy, by presenting FPMD results for the vis- 
cosity of liquid aluminum near the triple point. We then 
illustrate the potential importance of this type of calcu- 
lation by showing how FPMD calculations can be used to 
calculate the viscosity of liquid iron in the Earth's outer 
core. This is one of the key quantities in the theory of the 
Earth's deep interior, but also one of the most uncertain. 

There are basically two ways of calculating transport 
coefficients by simulation. The first is to subject the sim- 
ulated system to an explicit external perturbation (e.g. 
a shear flow or a temperature gradient) and calculate the 
steady-state response. Alternatively, one can apply the 
Green-Kubo (GK) relations, i.e. the relations between 
transport coefficients and correlation functions involving 
fluxes of conserved quantities ||]. The shear viscosity 77, 
for example, is given by: 



V 



(1) 



where (•) denotes the thermal average, V is the volume 
of the system, T is the temperature, fcs is the Boltzmann 
factor, and Pxy is the off-diagonal component of the stress 
tensor (a and /3 are Cartesian components) . The rel- 
ative merits of the two approaches have been much dis- 
cussed, but the GK method has the virtues of simplicity 
and ease of application, and will be used here. 

For a simulation having a given duration, single parti- 
cle properties, like the diffusion coefficient, can be calcu- 
lated more accurately than collective properties, like the 



viscosity. In the former case the statistical average can be 
done over time and over the particles, while in the latter 
one loses the possibility of averaging over the particles. 
To obtain the same statistical accuracy, collective proper- 
ties need much longer runs than single particle properties 
by a factor proportional to the size of the system. 

The first principles calculations presented here are 
based on density functional theory, pseudopotentials and 
plane-wave basis sets. The electron-ion interaction is de- 
scribed by ultrasoft Vanderbilt pseudopotentials (PP) [^j. 
Non-linear core corrections Q arc included throughout 
this work. The calculations have been performed using 
VASP (Vienna ab initio simulation package) |^ . In VASP 
the electronic ground state is calculated exactly (within 
a self-consistency threshold) at each MD step, using an 
efficient iterative matrix diagonalization scheme and a 
Pulay mixer Within this approach to FPMD it is 
important to provide a good starting electronic charge 
density at each time step, so as to reduce the number 
of iterations to achieve self-consistency. We have used 
the following charge extrapolation scheme: at the be- 
ginning of each time step the electronic charge density 
is extrapolated using the atomic charge density and a 
quadratic extrapolation on the difference, i.e. the charge 
is written as: p{t) — Pat{t) + 6p{t) where p{t) is the 
self-consistent charge density at time t, and Pat{t) is the 
sum of the atomic charges. At time t + dt the charge is 
written as the sum of the atomic charges, which can be 
calculated exactly and cheaply, and a quadratic extrapo- 
lation on 6 p. This scheme provides a much better start- 
ing charge than the standard quadratic extrapolation of 
the whole charge |Q. The electronic levels are occupied 
according to Fermi statistics corresponding to the tem- 
perature of the simulation. This prescription also avoids 
problems with level crossing during the self-consistent cy- 
cles. Forces and stress tensor are calculated using the 
Hellmann-Feynman theorem. The temperature is con- 
trolled using a Nose thermostat 

Our calculation of the viscosity of liquid aluminum was 
performed using the local density approximation (LDA) 
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for the exchange-correlation energy and a plane-wave 
cutoff of 130 eV. The temperature was 1000 K (70 K 
above the melting point) and the pressure was close to 
zero. We used a cell containing 64 atoms with periodic 
boundary conditions and T point sampling only. To inte- 
grate Newton's equation of motion for the ions we used 
the Verlet ||l^ algorithm with a time step of 3 fs, and the 
convergence threshold on the total energy was 1.5 x 10~^ 
eV/atom; with these prescriptions the drift in the to- 
tal energy was « 0.2 meV/atom per ps. The length of 
the simulation was 80 ps. The simulation was split into 
three parts, of lengths « 40 ps, « 25 ps, and w 15 ps, 
and at the beginning of each part the Nose thermostat 
was reinitialized, so as to avoid possible problems with 
the accumulation of the drift in the total energy. The cell 
volume was adjusted to give zero pressure; this resulted 
in a density oi p = 2470 kg m"'^, which is about 5% 
larger than the experimental one [ pT| . This discrepancy 
is probably due mostly to the LDA. 

There are five independent components of the traceless 
stress tensor, Pxyi Pyz^ -Pzx: 2{-^xx ^yy)^ ^zz)i 
the last two being equivalent to the first three by rota- 
tional invariance, so one can construct five independent 
stress autocorrelation functions (SACF). These are cal- 
culated by averaging Paf}{tA-tQ)Pai3{to) over time origins 
io- Since reinitialization of the thermostat causes small 
discontinuities, this averaging is restricted so that Iq and 
t + tQ never span a discontinuity. 

In Fig. |l] we display the average of the five SACF's 
divided by its value at t = for the aluminum system. 
Since the traceless part of the stress tensor has zero av- 
erage, goes to zero as t — > cx). The statistical error 
on (f){t) for all values of t is ss 1.5% of the value dX t — 
and after 0.3 — 0.4 ps the magnitude of 4){t) falls below 
that error. 

In Fig. ^ we display the integral dt'(l){t') of as 
a function of time. The limit value of the integral for 
i ^ cxD is the shear viscosity. The error that one makes 
in evaluating that integral grows with time, since one 
integrates the noise together with 4){t). We estimated the 
error in the integral as a function of time using the scatter 
of the SACF's calculated by splitting the simulation into 
many short intervals. We combined this estimate with 
an analytic expression for the error, and the resulting 
error estimate is displayed in Fig. ^. From the point 
where 4){t) falls below the noise one integrates only the 
latter, so one gains nothing by evaluating the integral 
beyond that point. If we assume that (j){t) is zero above 
t « 0.4 ps, we obtain the value r] = 2.2 ±0.1 mPa s, 
which should be compared with the experimental value of 
1.25 mPa s [Q. The fact that the statistical errors have 
been reduced to w 5% demonstrates that the viscosity of 
liquids can be calculated completely from first principles 
with quite high precision. Our calculated value of the 
viscosity differs significantly from the experimental one. 
In order to find a possible reason for this discrepancy 
we have simulated for 20 ps the aluminum system at the 
experimental density, p — 2350 kg m^^ (which gives a 



calculated negative pressure of 2 GPa). The results for 
the viscosity of this simulation are also reported in Fig. 
||. We obtain r/ = 1.4 ± 0.15 mPa s, which is in excellent 
agreement with the experimental data. 

A method which has sometimes been used in the past 
to obtain an approximate estimate of the shear viscos- 
ity exploits its connection with the self-diffusion coeffi- 
cient D via the Stokes-Einstein relation: _Dr/ ksT /2'Ka, 
where a is an effective atomic diameter. This relation 
is exact for the Brownian motion of a macroscopic par- 
ticle of diameter a in a liquid of viscosity r], but it 
is only approximate when applied to atoms. However, 
if a is chosen to be the radius ri of the first peak in 
the radial distribution function, the relation usually pre- 
dicts rj to within 40%. We have calculated D from our 
simulated liquid aluminum in the standard way from 
the long time slope of the mean-square displacement 
{\Yi{toA-t)—Ti{tQ)\^)^&Dt, as t ^ oo, where is the 
position of atom i. The resulting value D sa 5.2 x 10~^ 
m s^ , combined with the value a — ri = 2.65 A, gives 
77 « 1.6 mPa s, which agrees quite well with the Green- 
Kubo value. For the system at the experimental density 
we calculate D « 6.8 x 10~^ m'^s"^, which is consistent 
with the value of the viscosity. 

Although the statistical errors on r] are small, one 
might ask whether a 64-atoms system is big enough to 
calculate 77 accurately. The effect of system size on calcu- 
lated viscosities has been extensively studied for models 
such as the hard-sphere and Lennard-Jones liquids — close 
packed systems which are similar to liquid aluminum. For 
example Schoen and Hoheisel |l^ examined Green-Kubo 
calculations of the Lennard-Jones liquid for system sizes 
from 32 to 2048 atoms. They showed that even 32 atoms 
give meaningful values for 77 and that for 64 atoms the 
error on rj at the triple point is w 10%. We note that one 
indicator of system size effects should be the difference 
between 77 values calculated with the stress components 
Pxy etc. and those calculated from \{Pyy — Pzz) etc., 
since cubic periodicity breaks the spherical equivalence 
of the two types of SACF. We searched for this differ- 
ence in our calculations, but it is not large enough to be 
clearly visible above the statistical noise. 

Liquid iron under Earth's core conditions provides a 
good example of the importance of being able to calculate 
transport coefficients by FPMD. The Earth's liquid outer 
core consists mainly of molten iron, with light impuri- 
ties, of which sulfur is a likely candidate iQ. The shear 
viscosity of this liquid is a crucial factor in understand- 
ing convective circulation in the core, which is intimately 
linked both to heat transport and to the generation of the 
Earth's magnetic field. Nevertheless it is one of the most 
uncertain quantities in geophysics, with estimates from 
different experimental and theoretical methods spanning 
no less than 12 orders of magnitude |l^. A firmly based 
calculation of the viscosity is therefore important. In our 
recent FPMD calculations for liquid Fe ||l^,|l^ and liq- 
uid Fe/S 1^] under Earth's core conditions we used the 
Stokes-Einstein relation to obtain the estimate 77 = 13 
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mPa s. But the use of this relation under such extreme 
conditions can clearly be challenged. 

We present here the results that we have obtained from 
a direct calculation of the viscosity of an Fe/S liquid al- 
loy under Earth's core conditions. The calculations were 
performed at a thermodynamic state corresponding to 
the boundary between the solid inner core and the liquid 
outer core. In fact, the temperature at this boundary is 
uncertain, with estimates ranging from 4000 to 8000 K, 
while the density is quite accurately known to be about 
12000 kg m-3 The temperature in our simulation 

was set to 6000 K and the density was 12330 kg m^'^. 
The Fe/S mixture was taken to have a sulfur mole frac- 
tion of 0.1875, in line with the maximum estimates for 
sulfur abundance in the core [jl^ . The simulated system 
contained 64 atoms (52 Fe and 12 S) in the repeated cell, 
with F point sampling onlv We have used the generalized 
gradient approximation pU for the exchange-correlation 
energy, a plane wave cutoff of 350 eV, a time step of 1 fs 
and a self-consistency threshold of 1.5 x 10^^ eV/atom. 
The total length of the simulation was 10 ps and the 
drift on the total energy was « 8 meV/atom per ps. A 
complete description of the results of this simulation is 
reported elsewhere |^]. 

In Fig. ^ we display the average SACF (j){t) for this 
system. This is noisier than for the aluminum system, 
the statistical error at any time being « 4% of 4>{t = 0), 
but it is still possible to produce a useful estimate for 
the viscosity. The (/)(t) has gone to zero after t ~ 0.2 ps 
(Fig. 3), and the time integral of (j>{t) up to this point 
gives 77 = 9 ± 2 mPa s (Fig. 4), which confirms the ap- 
proximate correctness of the value estimated using the 
Stokes-Einstein relation. The error has been estimated 
as in the aluminum case. This value of 77 is about 10 times 
larger than the viscosity of typical liquid metals at am- 
bient pressure. This result is important since it provides 
concrete support for the approximation commonly made 
in magnetohydrodynamic models that the outer core is 
an inviscid fluid ||2l|,^^ undergoing small-circulation tur- 
bulent convection [^3t rather than large-scale global cir- 
culation. 

We have chosen the case of liquid Fe/S to illustrate the 
geophysical importance of first-principles calculations of 
viscosity, but there are many other areas where calcu- 
lated viscosities would be valuable. The viscosity of the 
highly compressed hydrogen/helium mixtures in the deep 
interior of the giant planets is a case in point. There have 
been recent attempts to deduce this viscosity from mea- 
sured gravitational anomalies for Jupiter, and input from 
first-principles calculations obtained using the methods 
we have described would be very useful. Insight into the 
percolation of molten minerals in the Earth and the other 
terrestrial planets could also be gained by first principles 
methods. In a completely different context, the viscous 
flow of materials is crucial in a wide variety of industrial 
processes. 

Although we have focused here on the shear viscos- 
ity, first-principles calculations should also be feasible for 



other transport coefficients like the bulk viscosity, the 
thermal conductivity and chemical interdiffusion coeffi- 
cients, all of which contribute to the attenuation of sound 
in fluids. 

In conclusion, we have demonstrated the feasibility of 
calculating transport coefficients such as viscosity from 
first principles, and we have presented results for the vis- 
cosity of liquid aluminum which agree well with experi- 
ment. We illustrated the importance of such calculations 
by presenting results for the viscosity of liquid Fe/S under 
Earth's core conditions. These results support recent ap- 
proximate estimates which indicate that the viscosity of 
the Earth's outer core is in the region of 10 mPa s. This is 
far lower than the values sometimes inferred from seismic 
and other measurements, and has important implications 
for the understanding of circulation in the core. Finally, 
we have pointed out that the techniques employed should 
also be applicable to other transport coefficients. 
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FIG. 1. Average over the five independent components of 
the autocorrelation function of the traceless stress tensor (j){t) 
calculated for liquid Al at 1000 K at the density of 2470 kg 
m~^. Values of (j>{t) are normalized by dividing by <?i(0). 
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FIG. 2. Viscosity integral of the average stress autocorre- 
lation function and its statistical error as a function of time 
for liquid Al at 1000 K. Results are shown for the calculated 
zero pressure density p = 2470 kg m~^ and the experimental 
density p = 2350 kg m~^. 
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FIG. 3. Average over the five independent components of 
the autocorrelation function of the traceless stress tensor (j){t) 
calculated for liquid Fe/S under Earth's core conditions. Val- 
ues of 4>{t) are normalized by dividing by 0(0). 
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FIG. 4. Viscosity integral of the average stress autocorre- 
lation function and its statistical error as a function of time 
for liquid Fe/S under Earth's core conditions. 
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